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ABSTRACT 

We study the dynamics of a system of two super-Earths embedded in a protoplanetary 
disc. We build a simple model of an irradiated viscous disc and use analytical prescriptions 
for the planet-disc interactions which lead to migration. We show that depending on the disc 
parameters, planets’ masses and their positions in the disc, the migration of each planet can 
be inward or outward and the migration of a two-planet system can be convergent (which may 
lead to formation of a first order mean motion resonance, MMR) or divergent (a system moves 
away from MMR). We performed 3500 simulations of the migration of two-planet systems 
with various masses and initial orbits. Almost all of them end up as resonant configurations, 
although the period ratios may be very distant from the nominal values of a given MMR. We 
found that almost all the systems resulting from the migration are periodic configurations. 
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1 INTRODUCTION 


It is widely accepted that planets form around young stars in proto¬ 
planetary discs at distances larger than the ones observed in evolved 
systems. Gravitational interactions between a planet and a disc ex¬ 
cites density waves in the disc. The existence of the waves results 
in a non-radial component of the force acting on the planet and, as 
a consequence, the planet migrates inwards ( [Goldreich & Tremaine | 
|1979[|T980| ). Since these pioneering works, significant progress has 
been made in the theory of the planet-disc interactions (e.g., |Tanaka| 
|et al.|2002||Paardekooper et al.|2010[|2011| >, showing that the di¬ 
rection as well as the rate of the migration depend on the planets’ 
masses and the disc properties in a complex way. 

The smooth convergent migration of the planets together with 
the eccentricity damping naturally leads to systems with P 2 /P 1 ~ 
(p + \)/p ([Lee & Peale|2002 Snellgrove et al.|200l Papaloizou 


|& Sz uszkiewicz 2005}. Therefore, one could expect that the obser¬ 

vational sample of multi-planet systems is dominated by resonant 
configurations. Nevertheless, |Fabrycky et al.] ( |2014} showed that a 
histogram of period ratios of systems discovered by the Kepler mis¬ 
sion does not reveal clear maxima at the positions of the MMRs. 

There are several explanations of those features presented in 
the literature. One of the proposed mechanisms are the stochastic 
forces acting on the planets. The forces might origin from the tur¬ 
bulence in a disc. |Nelson| |2005j ) showed that the turbulent density 
fluctuations can be larger than the spiral wakes excited by the plan¬ 
ets of a few Earth masses. The stochastic forces might also result 
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from the interactions between the planets and planetesimals (e.g., 
|Chatterjee & Ford|2014| . The interactions can disrupt a MMR if 
a planetesimal disc has a mass of > 0.2 x the mass of the plan¬ 
ets. Stochastic forces can be also incorporated in a model heuristi- 
cally (e.g., [Rein|2012| ). In such an approach, one can parametrize 
the migration, circularization and stochastic force. For certain val¬ 
ues of the parameters, it is possible to reconstruct the histogram of 
the period ratios. Nevertheless, |Hands et al.| ( [2bl4|) showed that in¬ 
creasing magnitude of the turbulence does not have to reduce sig¬ 
nificantly the number of resonant systems with respect to all the 
systems, but it rather reduces the survival rate of the systems. 


Dissipation due to the planet-star tidal interaction is another 
mechanism which could be responsible for moving the systems 


out of resonances dPapaloizou & Terquem 

2010 

Papaloizou 2011 

Batygin & Morbidelli 2013, Delisle et al. 

2014 

Delisle & Laskar 

2014). The tidal dissipation of the mechanical energy causes the 


inward migration of a planet, which is (in a two-planet system) 
faster for the inner planet. Nevertheless, the tidal dissipation is fast 
enough only for short-period planets (P < 10 days, this value de¬ 


pends on the type of planets, [Lee et al.||2013| ). The divergent mi¬ 
gration can also result from the interaction between a planet and a 
wake produced by another planet in the same system ( |Podlewska-| 
Gaea et al. 2012, Baruteau & Papaloizou 2013). This mechanism is 
efficient if at least one of the planets opens a partial gap. 

|Goldreich & Schlichting| < [20T4} show that the capture of plan¬ 
ets into a MMR is permanent if the equilibrium eccentricity e eq < 
(m/M ) 1 / 3 , where m and M are masses of the planet and the star, 
respectively. The value of e eq results from the balance between ex¬ 
citation of e due to passing through MMR and damping the eccen- 
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tricity due to planet-disc interaction. If the damping is not strong 
enough, e eq can be too high and the capture into MMR is only tem¬ 
porary. In such a case resonant configurations should be rare. On 
the other hand, as shown by [Deck & Batygin] j2015| ) although a 
given system resides in a MMR only some limited time, shortly 
after leaving one MMR the system goes into another MMR. 

In this work we present another scenario of shifting systems 
away from MMRs. We show that in a standard 1+1D a—disc model 
with a realistic opacity law, the disc can be divided into several re¬ 
gions of convergent and divergent migration. During the evolution 
of the disc, these regions move inwards at the viscous time-scale 
x vis . The planets migrate in the disc at rates which are in general 
different from x vis . Therefore, during the migration the period ratio 
can decrease (the system may be trapped into MMR) or increase 
(the system moves out of MMR) depending on the current positions 
of the planets in the disc. The final state of the system (in particu¬ 
lar, the period ratio) depends on the masses and initial orbits of the 
planets. Variety of final period ratios can result from the migration. 
One of the possibilities are systems with P 2 /P 1 ~ (p-j- l)/p for 
p = 1,2,3 or 4. Systems with the period ratios between the nomi¬ 
nal values of two neighbouring first-order MMRs (like 2:1 and 3:2) 
are also possible. The migration can also lead to formation of hi¬ 
erarchical systems of period ratios >10. The common feature of 
almost all these systems is that they evolve periodically and at least 
one of the two resonant angles oscillates. Moreover, almost all sys¬ 
tems with P 2 /P 1 <2.12 have two resonant angles which librate. 

In section 2 we introduce a disc model. In section 3 we present 
analytical prescriptions of the planet-disc interactions and study the 
migration rate as a function of planet’s mass and the size of its or¬ 
bit. In section 4 we discuss the migration of a single planet. In the 
next section, the evolution of a two-planet system is studied. We 
show that sizes and positions of the regions of convergent and di¬ 
vergent migration depend on the planets’ masses and evolutionary 
state of the disc. In section 6 we discuss main features of the sample 
of 3500 systems. Section 7 is devoted to limitations of the model. 
Conclusions are given in Section 8. 


2 A MODEL OF A DISC 

We consider a geometrically thin, viscous, axially symmetric disc. 
We assume that the planets do not affect the disc evolution and that 
the vertical and radial structures of the disc are independent (1+1D 
methodology, [Garaud & Lin|2007| ). Below we describe the model. 


2.1 The diffusion equation 


In order to study the migration of planets down to very small dis¬ 
tances from the star, we use the angular velocity profile different 
from usually used Keplerian profile £2 K = (GM/r 3 ) 1 / 2 , where G is 
the gravitational constant and r is a distance from the star. In the 
so called boundary layer ( |Pringle|1981) , Q decreases rapidly from 
a value close to fl K down to the rotational angular velocity of the 
star £2* <C £2 k- The boundary layer width L is much smaller than 
the star’s radius R *. Assuming £2* = 0, we chose 


Q(r) — Q k 




( 1 ) 


This particular form of the angular velocity is convenient when one 
wants to apply it to the diffusion equation, which then reads: 


dZ 

di 


= 3(1-?) 


dr 2 


(vE) + 


9 

2r 




( 2 ) 


where £ = (4/3)(/?*/r)*/ 2 and v is the turbulent viscosity coeffi¬ 
cient averaged over the vertical axis z and t is time. Far from the 
star £ —)• 0 and Eq. § tends to usually used formula. 


2.1.1 Photoevaporation 


The right-hand side of Eq. needs to be completed with a pho¬ 
toevaporation term Z w . We follow Matsuyama et al. (2003), but 
change slightly their Z w in order to have smooth transition between 
region of the disc from which the material is removed to the region 
of gravitationally bounded gas. We have 

*(*r- 


r>r g 


Z w — 


Eq exp 


(9?)' 


(1 -P)r g <r<r g . 


(3) 


where So = 1.257 x 3 ^ 2 gcm _2 s _1 , 71,4 is the 

temperature of ionized HII atmosphere above the disc divided by 
10 4 K, <l> 4 i is the ionizing photons luminosity divided by 10 41 s _1 , 
mo = M*/M 0 and r g = 1.26 T~l AU, . We adopted <T> 4 i = 1, 
(3 = 0.2 and 7^4 = 1. Equation ^ together with Z w is being solved 
with the following boundary condition. At r = r inner we assume con¬ 
stant accretion rate (Lv = const), while the outer boundary condi¬ 
tion reads Z(r outer ) = 0. 


2.2 The vertical structure equations 


To solve Eq. [2] one needs to know v at a given distance r and a 
given time t. This is achieved by solving standard equations of the 
vertical structure. The only modification we make is that £2(r) is 
given by Eq. |lj instead of £1 K - The equations read as follows: 


dP 

dz 

dF 

dz 

dT 

dz 


-pgz = -?®iz(^+ z -^j , 

£>(r) = pv(r^ =^pvn|(l-Q 2 , 


(4) 

(5) 

( 6 ) 


where P, p,F and T denote pressure, density, energy flux and tem¬ 
perature, respectively. The viscosity coefficient is given by the a 
prescription ( [Shakura & Sunyaev|1973[|D’Alessio et al.|1998) 


V=v(z) 


a P(z) 
p(z)£l’ 


a — const. 


(7) 


We assume that the energy is transported by radiation, convec¬ 
tion and turbulence ( jD’Alessio et al.||l998| >. The gradient V = 
d InT/d InP is computed with the help of the mixing length the¬ 
ory (we follow the approach of |Heinzeller et al.|2009| ). 


2.2.1 Boundary conditions 

Following ( [Papaloizou & Terquem||1999] > we write the boundary 
conditions, P s = P(zoo), F s = F(zoo), T s = T(zoo), F 0 = F(z = 0): 


Ps 

Ps 

0 

Po 


77 x ab 




0, 


8 k s jLim H 


-r, 


( 8 ) 

(9) 

( 10 ) 

(11) 
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Figure 1. Panel (a) Surface density profile for several values of constant M (labelled accordingly). Panel (b) Surface density evolution. Red curve is for the 
initial disc, black curves show E—profiles every 10 5 years. 


where Zoo is a height above which no energy is produced, the disc 
height H is defined such that p (H) = p(0) exp(—1), M is the accre¬ 
tion rate, x ab is the optical depth above the disc, K s is the opacity at 
z = Zoo, T b is the background temperature, ju is the mean molecular 
weight, k B means the Boltzmann constant and m H is the hydrogen 
atom mass. 


2.2.2 Stellar irradiation 


We use a simple approach to incorporate stellar irradiation of the 
disc, i.e., we modify temperature at the boundary T s by adding T irr 


in a forth power, i.e. T 4 := T s 4 + T 4 . Hueso & Guillot 2005), where 
T irr has a form of ( |Ruden & Pollack 


1991): 


Tin- — T± 


■_2_ (R*y 1 (R* 
3n\r + 2 r 


H\ fdlnH 


d\nr 


-1 


( 12 ) 


It is a common approach to assume d\nH/d\nr = 9/7 (e.g., Mor- 
dasini et al. 2015). This assumption seems to be a reasonable 
choice, because it is an equilibrium solution for a disc in regions 
dominated by the irradiation ( [Chiang & Goldreich|1997| and, on 
the other hand, T irr <C T s in regions dominated by the viscous heat¬ 
ing. Nevertheless, the approximation does not have to be good in 
regions where the viscous heating and the irradiation are of similar 
magnitude. We resign from this simplification and write T irr as: 


Tirr — T* 


3^[yJ + 2 Sa [t 


2i 1 


(13) 


where h = H/r is the aspect ratio and = rdh/dr is an unknown 
to be found in a self-consistent way. 


2.2.3 Opacity 

We use opacity tables from Semen ov et al.| ( 2003]> instead o f usu¬ 
ally used formulae ( |Ruden & Pollack|1991| ). |Semenov et ah| ( |2003| ) 
computed K for variety of dust models (several models of the grain 
structure and topology as well as the chemical composition). We 
checked that the differences does not affect the disc profile sig¬ 
nificantly. In further study we use a model of composite spherical 
grains of typical mass fraction of metallic iron. 


2.3 Solving the equations of the vertical structure 

We solve Eqs. {4])-([6]) for a given r and M starting from z = Zoo 
down to z = 0. Our unknowns are Zoo and s^. A method to solve 
the equations is the following. We choose some value of s^, let’s 
call it ^ adopted) . Once we have at a given radius r for a given M, 
we know T s . Then we find such Zoo for which Fq = 0. Then all the 
physical quantities (P, F , P, p) are known as functions of z, which 
gives us also H. Repeating the process for radii r — 5r and r + 5r 
we find = ^ compulcd) . A proper solution of the vertical structure 
equations means ^ computed) = ^(adopted) j n p r j nc jpi e? there are three 
possible types situations for a given M and r, i.e., an irradiated 
disc, a shadowed disc and a non-unique solution. Figure [T^ shows 
E-profiles for several values of M ^ lO -7 M 0 yr -1 . For a refer¬ 
ence, the disc of M = lO _9 M 0 yr _1 (red curve) is shadowed for 
log 10 r[AU]e[~0,~0.3]. 


2.4 Time evolution of the disc 


To solve Eq. {2} together with the photoevaporation term, Eq. {3}, 
we need to know v at given t and r. Nevertheless, it is not neces¬ 
sary to solve Eqs. 0-© at each time-step of the evolution of the 
disc. For a given r, v is a function of M and because E is a mono- 
tonically increasing function of M at a given radius, one can write 
v=v(E,r). This way we close Eq. j2| in a sense, that our only un¬ 
known is E = E(t, r) and v is given by E at given r and t. One needs 
to solve the vertical structure equations once for different r and M 
values taken from a grid (log 10 r[AU] G [—2,2] with a step of 0.005, 
log lo M[M 0 yr _1 ] G [—16, —6] with a step of 0.1). Exact values of 
v for given E are obtained with the help of the interpolation. 

The disc model is computed for the following values of the 
parameters: a = 0.004, M* = 1M 0 , R * = 2R 0 , the effective tem¬ 
perature of the star T e ff = 4000K. The initial profile Z(t = 0, r) oc 
(r/a 0 )~ b exp[-(r/a c ) (2_&) ] is taken after [Fortier et al. (2013]), 


where we chose b = 0.8, a^ = 5.0AU,ac = 10.0AU. The initial 
mass of the disc is 0.04 M 0 , the inner radius 0.02 AU, and the outer 
radius 100 AU. The outer radius as well as values of Z?, ao and ap are 
not crucial for the evolution of the inner parts of the disc, r < 10 AU. 
We integrate Eq. {2} until the disc is dispersed. Figure[lJ> shows the 
evolution of E. Red curve is for the initial profile. Black curves 
show the profiles every 10 5 yr. The shadowed disc region exists in 
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Figure 2. The migration rate as a functions of m and r. Green colour denotes the inward migration, while red is for the outward migration. The darker the 
colour is, the faster is the migration. The accretion rate and the eccentricity are given at each panel. 


first ~ 2.3 Myr of the evolution. The photoevaporation is the most 
effective for radii of a few AU. Similar evolution of Z-profile (apart 
from the shadowed region) can be found in ( [Mordasini et al.|20T2| . 


3 PLANET-DISC INTERACTION 

Because all the studies leading to analytical prescriptions of the 
torque are devoted to single planet-disc interaction (the density 
waves in a disc excited by one planet, affect only this planet, not 
the other planets embedded in the same disc) and finding the for¬ 
mulae for a general case is beyond the scope of this paper, we use 
here this approximation. We also assume that the planet does not 
affect the disc evolution. It means that once we solved Eq. [2] we 
use the solution L(t,r), T(t,r), etc. to compute the force acting on 
a planet. The disc structure is given on a uniform grid of log 10 r and 
with a constant output time-step of 5 kyr. Values of Z and T on the 
grid give local power indices yi and 72 (Z oc r _Yl , T oc r - ^ 2 ). The 
physical quantities at a given r and t are obtained with the help of 
the spline interpolation in log 10 r and the linear interpolation in t. 

The total torque acting on a planet moving in a circular orbit 
is a sum of the Lindblad torque T L and the corotation torque r c , 
T = T l + r c . We use formulae from ( [Paardekooper et ah|201lj ), 
where they study non-isothermal Type I migration. T is a function 


of masses of the star and the planet as well as the disc parameters 
(see |Paardekooper et al.|2011] for details). 

Once we compute T for a planet of a given mass ra in an orbit 
of a given size r, the acceleration acting on the planet is given by 


_ r _ fhy/jLia _ L 0 

Jt ~ 2x a ’ Za ~ 2F 2F 


(14) 


where \ a is the time-scale of migration, a is the semi-major axis, 
ju = G(M * -bra), ra = (1/M* + 1/ra ) -1 is the reduced mass of the 
planet, Lq is the angular momentum of the planet in circular orbit 
and r is the astrocentric velocity of the planet. 

Because in general the orbit is non-circular, we take into 
account the eccentricity waves (Tanaka & Ward 2004). The for¬ 
mulae in that paper are given in a specific reference frame. The 
planet starts from the pericenter and the size of the orbit is fixed. 
The force acting on the planet is given as an explicit function 
of time. The formulae cannot be used directly in a more general 
case. Assuming coplanar orbits, we rewrite their Eqs. [38]—[40] as 
F r =A (0.057x + 0.176y) and F§ =A (—0.868x + 0.325y), where 
A = mM^ 2 C~ 4 a 6 eLQ. 6 , e is the eccentricity and C s is the isother¬ 
mal sound speed. Cartesian coordinates of the planet (v,y) are given 
in the orbital reference frame (the x-axis points the pericenter). The 
acceleration acting on the planet reads: 


fe = F*ir+ 



r 

r 5 


L=\\rx r II 


(15) 
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When the orbit is eccentric, the corotation torque is weaker 
than for a circular case ( [Fendyke & Nelson|2014| ). The corotation 
torque for the orbit of eccentricity e is given by the torque for a 
circular orbit T c (0) multiplied by a factor of a value smaller than 
unity, T c (e) = exp(— e/ef) r c (0), where ef = h /2 + 0.01. 

For masses of the planets > 1OM 0 and for not very massive 
disc, L(1AU) < lOOgcm -2 , formulae of Type I migration may not 
describe the disc-planet interaction correctly. Therefore, we use a 
simple analytic prescriptions from ( [Dittkrist et al.|[20l4] ) for the 
transition between Type I and Type II migration formulae. 

The last term in our model is the gravitational force originating 
from an unperturbed disc. Papers devoted to the planet-disc inter¬ 
actions are usually focused on a and e evolution, thus the axially 
symmetric component of the potential from the disc is not taken 
into account as it leads only to the rotation of pericenter. When 
there are more planets in the disc, changes of the frequencies of 
the apsidal lines rotation affects the dynamics of the system. When 
Z(r)-profile is given at time t , the central force acting on a planet 
has a form of f r = — d(r)r, where d(r) stems from integration of 
the planet-disc interaction over a whole axially symmetric disc. 

Formulae for F r , and T c (e) depend explicitly on e and a. 
The quantities have a meaning of the osculating Keplerian ele¬ 
ments. Nevertheless, when the disc is massive and the planet is 
moving in a relatively wide orbit, the axially-symmetric component 
of the force from the disc can be a significant fraction of the central 
force from the star. A planet in a circular orbit, i.e., r(t) = const, 
is moving with a super-Keplerian velocity. Without taking into ac¬ 
count the additional force from the disc, one obtains e > 0 (also 
a will be incorrect) when uses the standard formulae for the os¬ 
culating Keplerian elements. The value of e one obtains can be 
even >0.1, which leads to incorrect values of f e and T c . There¬ 
fore, when computing a and e we take // = ju + d(r) r 3 instead of 

Id. 

Finally, the equations of motion for i -th planet read as follows: 
n = -mn +FI +f T +f e where 


/('') = - Y 

J p-p 


Grrii 


Gmi 


1 




r 3 r J 
fri r J 


(16) 


is the A-body component. The formulae given in this section are 
being used to integrate the equations of motion of the system of 
two planets. Before we go to this part of the analysis, in the next 
section we study how the migration rate depends on the planet’s 
mass and position in the disc of given constant accretion rate. 


3.1 The migration rate as a function of a and m 

Figure[2^,b,c present the migration rate for a planet in a circular or¬ 
bit embedded in a disc of constant accretion rate of 10 -8 ,10 -9 and 
10“ 10 M 0 yr _1 , respectively. The migration rate is higher for larger 
M and lower for smaller M. For M = 10 -8 and lO _9 M 0 yr _1 , 
there exist regions of the outward migration in the shadowed part 
of the disc (they are limited to mass ranges between ~ 1 and 
~ 1OM 0 ) as well as close to the star (this regions extends for a 
whole range of planet’s masses). For M = 10 -9 M Q yr -1 the out¬ 
ward migration occurs also for r < 0.1AU. For lower accretion rate 
M = lO -lo M 0 yr -1 almost all the (a, m )-plane corresponds to the 
inward migration. The migration rate depends also on e. Figure |2fl 
presents T^-map for M = lO -9 M 0 yr -1 and e = 0.02. All regions 
of the outward migration disappears (compare with Fig.[2|3). 

The existence of the regions of outward migration is important 
for planets growth during the migration (e.g., |Bitsch et al.|2014| ). 



Figure 3. Colour-coded: Migration rate of a single planet (m = 7 M®, e = 0) 
embedded in the disc as a function of time and the distance from the star. 
Vectors: Velocity field of a planet. Blue curves show the evolution of the 
systems of chosen initial semi-major axes (white symbols at t = 0.5 Myr). 


For a single-planet system, the eccentricity is damped to zero and 
this region exists. The result is that the planet is halted at the dis¬ 
tance of a few AU until the mass grows enough to omit the trap. 
However, when two planets migrating in a disc form a mean mo¬ 
tion resonance, their eccentricities are excited and the planets can 
omit the barrier even if their masses are small. This might lead to 
qualitative difference in the evolution of orbits as well as the mass 
growth of a single-planet system with respect to a multi-planet sys¬ 
tem. We leave this problem to future works. 


4 SINGLE PLANET MIGRATION 

In the previous section we discussed the migration time-scale of 
a planet of given mass in a disc of a given accretion rate. A se¬ 
ries of migration maps obtained for decreasing M can be related 
to different evolution stages of the disc. The disc evolves approx¬ 
imately as a M = const-disc, apart from the late stage of evolu¬ 
tion, when the photoevaporation plays an important role. Figure [3] 
presents x a —maps at (f,log 10 fli)— plane for m = 7M 0 and e = 0. 
The outward migration region occupies relatively large range of 
radii and exists up to ~ 2.3 Myr of the disc life-time. The blue 
curves show the evolution of single-planet systems for different ini¬ 
tial orbits. A planet which starts the evolution in or above the red 
region of the map cannot reach orbits of sizes significantly smaller 
than 1 AU. The planet reaches the trap (T = 0, % a °o) after a short 
time and migrates further inward at the time-scale of the disc evolu¬ 
tion, which is ~ 1.6 Myr. On the other hand, when a planet starts be¬ 
low the red area, it can migrate inward faster than the disc evolves. 
A planet starting at r ~ 0.8 AU can reach r ~ 0.02 AU. 

Before going to the next section, in which the evolution of 
two-planet systems is studied in more details, one can make an 
easy observation from the inspection of Fig. [3] When two plan¬ 
ets are initially located outside the trap and are migrating conver- 
gently, they can form a MMR. The convergent migration means 
that before forming the MMR, the outer planet migrates faster than 
the inner one (Ti > 12 ). When the planets migrate as a resonant 
pair, the outer planet migrates slower and the inner planet migrates 
faster with respect to the non-resonant case. It means that there 
has to be some torque exchange between the planets. The torque 
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Figure 4. Colours: Time-scale of the period ratio variation, Zx-Arrows: Vector field, i.e., the x-component of a vector is given by X multiplied by a time interval 
of twee = 1.5 x 10 4 yr, the y-component is given by d(log 10 a\)/dt multiplied by f V ec- Green colour means that Zx > 0 (convergent migration), red colour denotes 
Zx < 0 (divergent migration). The darker the colour is, the larger is \zx\- Masses of the planets are m\ = m 2 = 7M® (panel a), m\ = 6 M®,m 2 = 8 M® (b) and 
m\ = 8 M® ,m 2 = 6 M® (c). For all these plots t = 5 x 10 5 yr. Panel (d) shows the view of x^-map for m\ = 8 M® ,m 2 = 6 M® and t = 1.5 x 10 6 yr. 


acting on the inner planet, rp res , decreases, and the one for the 
outer planet, l^res* increases when compared to the non-resonant 
case, i.e., r^ res <C Ti, non-res an d T2,res ^ ^2,non-res* The torque act¬ 
ing on a planet which resides in the red region of the map is pos¬ 
itive in absence of the resonant interactions with the outer planet, 
Fi,non-res > 0- It is possible that ri res < 0, thus the inner planet can 
migrate inward and pass through the trap. The same cannot happen 
for the outer planet, because if r2 5non -res > also I^res > 0- In such 
a situation, the outer planet stays at a trap and the system becomes 
hierarchical ( Pi/P\ 1). The trap plays a very important role in 
forming the architecture of two-planet systems. In configurations 
with more than two planets, a ’’convoy” of planets can lead to a sit¬ 
uation in which all the planets, apart from the outermost one, pass 
through the trap and reach inner regions of the disc. We will study 
systems with larger number of planets in the next paper. 

Another important feature of the T a — maps is that at a given 
epoch T a has minima and maxima along log 10 a. The divergent mi¬ 


gration occurs when < x a ^ (if both xs are positive). It means 
that there are orbits for which a system, which initially resides in 
MMR, starts to diverge out of the resonance. It happens when the 
inner planet is close to the local minimum of z a . Although such a 
simple reasoning is valid for a case of equal masses m\ = m 2 , it 
shows that regions of divergent migration should exist in the disc. 


5 EVOLUTION OF TWO-PLANET SYSTEMS 


If T a i and % a ,2 are the time-scales of the migration of the inner and 
the outer planets in a pair, respectively, the time-scale of the period 
ratio variation reads as follows: 


ix = - 


X 

3P 


v _ p 2 2 Ta : iZa : 2 

X = —, Tx = - -• 

p l d 


(17) 


If both x a? i and 1:^2 are positive and Vi > x «,2, than %x > 0 and 
the period ratio decreases. If x a ^\ < the migration is diver- 
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Figure 5. Evolution of a chosen initial system of m\ = 8M®, m 2 = 6M®, log 10 ai = 0.6, P 2 /P 1 = 2.6. Initial epoch t = 0 corresponds to t = 5 x 10 5 yr of the 
disc evolution. Plots of a and e evolution present the results for both planets, the inner one (red) and the outer one (blue). White area bordered with solid black 
curve in panel (a) denotes the transition region in a disc between shadowed and irradiated parts. This region corresponds to the outward migration. Remaining 
panels present the evolution of P 2 /P 1 and the angles defined at corresponding panel. Characteristic stages of the evolution of the system (different ranges of t) 
are marked with different colours (see the text for details). 


gent. Similarly to the migration maps presented in Fig. [3] one 
can make maps for Ty. Figures [4^, b and c present x^-maps on 
(P 2 /P 1 , log 10 a i ) —plane computed for t = 5 x 10 5 yr and three dif¬ 
ferent pairs of masses m\ = m 2 = 7M 0 (panel a), m\ = 6 M 0 , m 2 = 
8 M 0 (b) and m\ = 8 M 0 ,m 2 = 6 M 0 (c). If a two-planet system is 
located in a green region, will decrease, while if the system 

is located in red regions, the period ratio will increase. As colours 
give only the information on the P2/A evolution, the arrows are 
over-plotted on the maps to complete this information. If the arrows 
point down/up, the migration of the inner planet is inward/outward. 
Directions left/right correspond to convergent/divergent migration. 

The view of the (P 2 / P \, log 10 a i ) —plane depends on the plan¬ 
ets masses. When the inner planet is less massive than the outer 
one, larger fraction of the plane corresponds to the convergent mi¬ 
gration (see Fig.| 4 ] 3 ). When m\ > m 2 , most of the plane corresponds 
to the divergent migration (see Fig. |4]e). One might expect that for 
mi < m 2 MMRs form more frequently than for mi > m 2 . 

As the disc evolves, the view of the —maps changes. Fig|4jl 
presents the map for m\ = 8 M 0 , m2 = 6 M 0 at t = 1.5 x 10 6 yr. 
Also the position of the system in the plane changes (compare po¬ 
sitions of an example system plotted with white symbols at panels 
c and d). It means that during the evolution, the system may reside 


in a convergent/divergent region some limited time and the period 
ratio may increase or decrease depending on the evolution stage. 
The final configuration depends in a complex way one the initial 
orbits, planets’ masses and disc evolution time scale. 

5.1 An example of the evolution of a two-planet system 

Figure [5]presents the evolution of an example system. Both planets 
of masses mi = 8 M 0 and m 2 = 6M 0 start the migration outside the 
trap (see Fig. [5ji, the white region bordered with black solid curve 
denotes the shadowed/irradiated disc transition region, see also the 
position of the system plotted with white symbol in Fig.[4^) at the 
epoch 0.5 Myr of the disc evolution. Ranges of time, corresponding 
to characteristic behaviour of the system, are colour-coded. 

Initially the planets migrate convergently, the period ratio de¬ 
creases from 2.6 down to ~ 1.65 (Fig.[5^). The convergent migra¬ 
tion is very fast, thus the system crosses 2:1 MMR. It does not reach 
3:2 MMR, because before that it moves into the region of the diver¬ 
gent migration. The period ratio increases approaching 2. The first 
critical angle of 2:1 MMR = 'ki — 2 X 2 + G5i (Fig. [ 5 ) 3 ) librates 
when P 2 / P\ is sufficiently close to 2. This period of the evolution 
is coloured with light red. At t ~ 0.9 Myr the inner planet reaches 
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Figure 6. Time-scale lx (colours) and a vector filed for the inner planet, 
which is involved in 2:1 MMR with the outer planet. m\ — m2 = 7 M®. 


the trap, while the outer planet keeps migrating inward. The result 
is that P 2 /P 1 decreases again and the system reaches 3:2 MMR. 
Before that, for t G (~ 0.9, ~ 1) Myr, when P2/P1 is still close to 
2, both resonant angles of 2:1 MMR librate (the second angle is 
defined as — ^l — 2 A ,2 + 052 , see Fig.J^Jl). This fragment of the 
evolution is coloured with light green. When both critical angles 
librate, also the difference between the longitudes of pericenter, 
A05 = 05i — 052, librates (Fig.[5j*). The period ratio keeps decreas¬ 
ing and when P 2 /^i ~ 1.7 the critical angles of 2:1 MMR start 
to rotate and the angles of 3:2 MMR ((^ = 2?ti — 3?L2 +®i and 

(j)^ = 2 A 4 — 3 X 2 + 052 , see Fig. |,h) start to librate. This fragment 
of t-axis, which corresponds to trapping the system into 3:2 MMR 
is coloured with light violet. At t ~ 1.1 Myr P 2 /P 1 approaches 1.5 
and the amplitudes of the librations of both resonant angles are very 
small. The eccentricities reach their equilibrium values <0.01. The 
planets migrate as a resonant pair for ~ 0.5 Myr (light blue colour). 
The inner planet passes through the outward migration region at 
t ~ 1.7 Myr. The outer planet reaches the trap and migrate at the 
viscous time-scale. The inner planet migrates faster, thus the sys¬ 
tem goes out of 3:2 MMR. The resonant angles of 3:2 MMR librate 
until P2/P1 < 1 7, after which the angles of 2:1 MMR start to librate 
again. When P2/P1 passes through 2 , the behaviour of the angles 
changes. The first angle librates around 0 instead of 180 degrees, 
and the second angle librates around 180 instead of 0 degrees. The 
second angle librates only when P2/P1 < 2 . 12 . 


fc (l) 


oscillates. Both 


The final period ratio is ~ 2.4 and only 
eccentricities are very small ~ 10 -4 . As shown in (Batygin & Mor- 


bidelli 2013), for very low e the frequency of the apsidal motion can 
be high enough to make the critical angle oscillate even for P2/P1 
significantly distant from (p + l)/p. Papaloizou ( |201 1) also stud¬ 
ied resonant behaviour of a system far from the nominal position 
of the resonance. As the system moves towards higher P 2 /A, its 
eccentricities are very small, but the resonant angles keep librating. 


5.2 Conditions for keeping a system in MMR 

Once a system is locked in MMR, one can ask at what radii and 
what evolution stage of the disc this configuration can be preserved. 
As we illustrated in the previous section, planets can remain in the 
region of the convergent migration only for some limited time, de¬ 


pending on the time-scales of their migration and the disc evolution. 
Figure [ 6 ] shows Tv-maps for 2:1 MMR for m\ = m 2 = 7M 0 . For 
a given a\, a value of 0,2 is such that P2/P1 = 2 . If initially reso¬ 
nant configuration is placed in the red region, the system will move 
out of the resonance. If the system is placed in the green region, its 
period ratio will be constant (although the green colour means that 
P2/P1 decreases, it holds true only for a non-resonant case). 

If the system is non-resonant, the evolution of the semi-major 
axes of the planets are given by formulae: 


^(non-res) 


2aiY\ 
£.1.0 ’ 


. (non-res) 
a 2 


2a 2 T 2 

^ 2,0 


(18) 


where F\ is the torque acting on the inner planet and +2 is the 
torque acting on the outer planet. The angular momenta of the plan¬ 
ets moving in circular orbits are given by L\ q and Z, 2 5 o for the inner 
and the outer planet, respectively. When the system is locked in 
MMR, the ratio of the semi-major axes is constant. Another con¬ 
strain is that the sum of the torques acting on the planets is the same 
in both resonant and non-resonant situations. Therefore we have the 
following formulae for a in the resonant case: 


.(res) _ 2«1 (ri +r 2 ) (res)_ 2 a 2 (ri+r 2 ) 

— J . J 1 a 2 — T I t • 

L \,0 + £2,0 £1,0+£2,0 

The vector field for a\ (shown in Fig. [ 6 } corresponds to the res¬ 
onant case. As we can see, the inner planet migrates faster than 
the disc evolves. Sooner or later an initially resonant configuration 
moves into the region of the divergent migration. At t ~ 2 Myr the 
orbits of log 10 fli[AU] ~ — 0.4, E (— 1 , —0.5),—1.5 correspond 
to the divergent migration, and systems located there cannot re¬ 
main in 2:1 MMR. At t > 3 Myr the resonance can be sustained 
only for r < 0.1AU. When the inner planet is less massive than the 
outer one most of the plane corresponds to the convergent migra¬ 
tion and the 2:1 MMR can be formed and preserved for almost all 
radii log 10 <zi[AU] < —0.5. When m\ > m 2 , most of the 'lx —map 
corresponds to the divergent migration. In such a case, a resonant 
system is unlikely to form due to migration. 


5.3 Migration for different masses and initial orbits 

As we showed above, the final orbital configuration depends on the 
masses of the planets as well as on the initial orbits. Figure[7]shows 
trajectories of systems starting from different initial positions at the 
(jyPi,log 10 ai)-plane. Panel (a) shows the results obtained for 
m\ = m 2 = 7M 0 . Initial positions of sixteen systems are marked 
with white symbols over-plotted at Tv—map computed for the ini¬ 
tial epoch. Eight of the systems (group 1) initially reside above the 
red region, another five (group 2 ) are located within this region, 
while three systems (group 3) below the region. The systems from 
the first group initially tends towards 2:1, 3:2 or 4:3 MMR. The 
planets migrate slightly faster than the disc evolves. They reach the 
red region, i.e., the inner planet passes through the outward migra¬ 
tion region and migrate further, while the outer planet is halted at 
the trap (see the discussion in Section 4). The systems start to di¬ 
verge from the resonances. After the trap vanishes at t ~ 2.3 Myrs 
(see Fig. [3j, the systems are again located in the convergent mi¬ 
gration zone. However, because the disc is of low mass and shortly 
after that moment it disperses, the systems do not have enough time 
to reach 2:1 MMR. The final P2/P\ are between 2 and 3. 

The systems from the second group evolve towards large 
P2/P1 • In this case, the inner planet in each system is initially below 
the trap and migrates inward, while the outer planet is stopped at 
the trap. Until the trap disappears letting the outer planet to migrate 


© 2012 RAS, MNRAS 000,|T][T2] 

























Migration of two planets in a disc and formation ofMMRs 9 



1 1.5 2 2.5 3 3.5 4 1 1.5 2 2.5 3 3.5 4 

P 2 /Pl P 2 /Pi 


Figure 7. The evolution of 16 initial systems (white dots) plotted with black solid curves. Final configurations are marked with green dots. Planets’ masses are 
given in top right-hand corner of each panel. The Zx ~maps for the beginning of the simulation is colour-coded. Vertical dashed lines denote MMRs. 


inward, the inner one is already very close to the star. These sys¬ 
tems end up as hierarchical configurations. Their final period ratios 
are larger than 4, thus they are not shown in the plot. 

The systems from the third group reach the resonances (2:1, 
3:2 or 4:3) and because they migrate faster than the disc evolves, 
they avoid the trap. However, at some moment of the evolution, the 
systems reach the second region of divergent migration (see Fig. [6] 
the red band initially located at a\ ~ 1AU and moving inward 
down to a\ ~ 0.1 AU). For some time (when a given system re¬ 
sides in the divergent migration region) their period ratios increase. 
For instance the system initially located at (P 2 /P 1 ,log 10 a\ [AU]) = 
(1.4,0) goes out of 4:3 MMR, and pass through 3:2 MMR. Af¬ 
ter the systems leave the region of divergent migration, they go 
towards MMRs (not necessarily the ones they left). After being 
locked in MMRs, the systems migrate inward as resonant pairs un¬ 
til they reach the third region of divergent migration (see Fig. [6] at 
log 10 ai [AU] ~ — 1.5 att > 2 Myr). They again diverge slightly out 
of the resonances and return. The final configurations are resonant 
and short-period (Pi in each pair ranges from one to a few days). 

Different behaviour can be observed in Fig.[7j). When mi > m 2 
most of the map is red (see Fig. El ,d) and none of the systems 
ends up as resonant or close-to-resonant configuration. For a refer¬ 
ence, the evolution of the system initially located at logi 0 a\ [AU] = 
0.6 and P 2 /P 1 = 2.6 was shown in Fig. [5] 


6 PERIODIC EVOLUTION 

We performed 3500 simulations for different planets’ masses and 
initial orbits. The masses of both planets are chosen from the 
range of [2,10] M®, initial log 10 ai[AU] G [—0.1,0.6] and initial 
P 2 /P 1 E [1.3,2.7]. Distributions of all these quantities are uniform. 
For almost all the final systems, at least one of the resonant angles 
of a given MMR oscillates/librates, and if we take into account only 
the systems with P 2 /P 1 < 2.12 (about half of the sample), for al¬ 
most all the solutions both angles librate. Only the first order res¬ 
onances are present. Even close to the nominal positions of higher 
order MMRs (5:3, 7:5, 8:5, etc.), the critical angles of these MMRs 
rotate. The evolution of a few representative systems which ended 
up in or close to 2:1 MMR are illustrated in Fig. [8] Each row of 
this figure is for one solution, they are enumerated with I, II, III 
and IV, respectively. The left-hand, middle and right-hand columns 
show the evolution of the semi-major axes, the period ratio and the 
resonant angles, respectively. 

The first example (I) is a system which migrates convergently 


down to 2:1 MMR and stays there during the whole evolution. 
The period ratio is slightly higher than 2, and the amplitudes of 
the resonant angles librations are very small. The first resonant an¬ 
gle of 2:1 MMR librates around 0, while the second one librates 
around 180 degrees. System II initially migrates convergently to¬ 
wards 3:2 MMR, stays close to the resonance for some time and 
then migrates out of this MMR. The final period ratio is ~ 1.75. At 
the end of the migration, both resonant angles of 2:1 MMR librate. 
The first one around 180 degrees, while the second one around 0. 
System III initially migrates towards 3:2 MMR, next it moves out of 
the resonance. It passes through 2:1 MMR, reaches P 2 /P 1 ~ 3 and 
returns. The final period ratio ~ 2.1. Both resonant angles librate, 
the first one around 0, while the second one around 180 degrees 
(like in system I). It is a generic feature of all the systems origi¬ 
nating from the migration within the model described in this paper. 
If a given system ends up with P 2 /P 1 below the nominal value of 
a given resonance, the first critical angle (which is always the one 
with 05 1 in its definition) librates around 180 degrees and the sec¬ 
ond one around 0. If P 2 /P 1 is larger than the nominal value, the 
situation is reversed. System IV ends up with P 2 /P 1 ~ 2.2 and only 
the first angle oscillates, the second one rotates. 

Almost all the solutions are characterized by oscila- 
tions/librations of at least one resonant angle. Another common 
feature of almost all the systems is that their evolution (after the 
disc is dispersed) is periodic. The final configurations whose migra¬ 
tion were illustrated in Fig.[8]were taken as initial conditions for the 
A-body simulations (see Fig. [9J. The integration time is 10 5 x P2. 
In all the cases, the evolution repeats every certain period. The pe¬ 
riod ~ 30 days for system I (it is almost exactly the period of the 
outer planet, as both the apsidal lines move very slowly, i.e., with 
a period of ~ 20 yr). The periods of the evolution of systems II, III 
and IV is shorter or longer than P 2 , depending on if P 2 /Pi is smaller 
or larger than 2. In all these cases the rotation of the apsidal lines is 
only a few times slower than the orbital motion of the outer planet. 
Eccentricities are also small ~ 10 -4 in these systems. Similar be¬ 
haviour of the angles can be observed for systems close to 3:2, 4:3 
and 5:4 MMR. 

As we observed, the resonant angles of a given first order reso¬ 
nance (see Fig. [8} oscillate/librate around 0 or 180 degrees even for 
systems with P 2 /P 1 significantly far from the nominal value of this 
MMR. It is possible because the apsidal motion for low eccentric 
systems may be very fast ( [Batygin & Morbidelli||2013] >. The ec¬ 
centricities are getting lower and the amplitudes are getting higher 
when P 2 /P 1 is moving away from a nominal value of a given MMR. 
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Figure 8. Migration of example initial configurations which ended up in or close to 2:1 MMR illustrated in semi-major axes-, period ratio- and resonant 
angles-plots (from left-hand to right-hand column, respectively). Each row is for one system (they are enumerated with I, II, III and IV). The masses of the 
planets for each configuration are given in the left-hand column panels for the reference. The critical angles of the resonance are (J)!, 1 } (red dots) and (j)^ (blue). 


Each system in the sample where integrated (within the A-body 
model) for 10 5 x P 2 taking as initial conditions the final orbital el¬ 
ements from the migration simulations. Figure [T0| shows how the 
amplitudes of librations of the resonant angles depend on the pe¬ 
riod ratio. The amplitudes for different angles are plotted with dif¬ 
ferent colours (see the caption to Fig. [Toj. We can see that there 
are values of P 2 IP 1 for which the amplitudes of two neighbouring 
MMRs tends to ~ 180 degrees. These values mean the borders be¬ 
tween two neighbouring first order MMRs (they are marked with 
vertical lines). For P 2 /P 1 € [~ 1.7, ~ 1.72] there are solutions for 
which and oscillate. This range of the period ratio can 
be considered as the border between 2:1 and 3:2 MMRs. Similar 
border can be seen between 3:2 and 4:3 MMRs (at P 2 /P 1 ~ 1.4). 
A border between 4:3 and 5:4 MMRs cannot be noticed because 
there is no solutions in this range of P 2 /P \. The outer border of 
2:1 MMR takes place at P 2 /P 1 ~ 2.12. The amplitude of the libra- 
( 2 ) 

tion of tends to ~ 180 degrees for this value of the period ratio. 

(2) 

In fact, for P 2 /P 1 >2.12 the amplitude for equals 360 degrees 
(it is not shown on the plot, as also the amplitudes for 3:2, 4:3 and 
5:4 MMRs equals 360 degrees in this range of the period ratio). An 
interesting observation one can make by looking at the amplitudes 
diagram on Fig. [TO] is that all the systems are located along well 
defined curves. It means that the amplitudes of the librations of the 
resonant angles do not depend on the planets’ masses. 






^ [deg] 4&i [deg] 

Figure 9. Evolution of the resonant angles of 2:1 MMR. Each initial con¬ 
dition were integrated for 10 5 x P 2 within the V-body model of motion. 
Labels I, II, III and IV corresponds to respective systems shown in Fig. [8] 
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Figure 10. Amplitudes of AG3 (black symbols) and the resonant angles of 
2:1 MMR - red, - light red), 3:2 MMR ((j)^ - green, (J)^ - light 

green), 4:3 MMR ((j)^ 1 ] - blue, (j)^ 2 ] - light blue) and 5:4 MMR ((j)^ 1 ] - violet, 

( 2 ) 

(|> 5-4 ~ light violet). Horizontal arrows show the ranges of first order MMRs. 

7 LIMITATIONS OF THE MODEL 

There are several important assumptions of the model, which we 
list below. We plan to improve the model in future works. 

7.1 Disc model 

A natural assumption of 1+1D model is that the radial and the az¬ 
imuthal structures of the disc can be studied separately. Moreover, 
we assumed that the turbulent viscosity is given by a prescription 
and that a is constant both in r and in z-direction. If the turbulence 
results from the magnetorotational instability (MRI |Balbus & Ha\w| 
|ley|1991[ |Hawley & Balbus|1991| >, the material of the disc needs 
to be sufficiently ionized to make this mechanism work efficiently. 
This condition is difficult to be fulfilled in cold dusty discs {bT aes 
|& Balbus||l994| . In such regions there might exist the so called 
’’dead zones” jGammie|1996| l, which do not transport the angu¬ 
lar momentum. The ionization can occur due to collisions between 
particles, which is efficient for T > 1000K and takes place near 
the star, r < 0.1 AU for relatively fast accretion rates of 10 -9 to 
10 -8 Mq yr -1 or even closer when the accretion is slower. Further 
from the star the galactic cosmic-ray ionization ( |Gammie|1996| ) can 
be effective in regions of low optical depth (low density regions or 
close to the disc surface). In the latter case the midplane regions are 
neutral or only slightly ionized. In more realistic model a should 
be a function of r and z (e.g., |Kretke & Lin|2010j ). 

7.2 Planet-disc interaction 

There are three important assumption we made about the interac¬ 
tion between the disc and the planets. We assumed that the disc 
evolves as it was not affected by the planets. For planets of a few 
Earth masses and relatively massive disc (M > 10 -9 M Q yr -1 ) the 
assumption is good. The formulae used in this paper takes into ac¬ 
count only the fact that a planet locally disturb the disc, i.e., it opens 
a partial gap (transition from Type I to Type II migration |Dittkrist| 
|et al.|2014| ). Second aspect is a single-planet approximation. One 
planet excites the waves in the disc and only this planets ’’feels” the 
gravitational force, which is a result of this perturbation. As shown 
in ( |Podlewska-Gaca et al.|2012[|Baruteau & Papaloizou|2013| > this 
effect may be important, at least for planets which are able to open 
a partial gap in the disc. Another simplification is that all formu¬ 


lae used here were obtained for models with given yi and 72 (the 
power indices of Z— and T -profiles), which were constant in suf¬ 
ficiently wide range of the disc (wide enough to encompass the 
wakes produced by the planet). Due to numerous transitions be¬ 
tween the opacity regimes, a range of r for which one can assume 
Yi and 72 to be constant may be narrower than the wake structure 
(in particular when the Lindblad torque is considered, as the coro¬ 
tation torque origins from a very tiny region around the position of 
the planet). On the other hand, the Lindblad torque is less sensitive 
to yi and 72 than the corotation torque is. 

7.3 Planet mass growth 

Another important issue is the planet mass growth. As discussed in 
( [Bitsch et al.|2014| ) when a planet reach a trap (which works only 
for planets in certain mass range) it can stay there for some time, 
accrete material from the disc and finally pass the trap. Moreover, 
even if there is no trap in the disc, the evolution of a planet whose 
mass is not constant is, in general, different from the evolution of 
constant mass planet (e.g., |Alibert et al.|2013[|Coleman & N elson 
|2014[ |Machida et al.||2010| . One could say that after some time 
the mass of the planet is almost constant. Still, one needs to start 
the integration from some initial orbits. On contrary, the model of 
growing mass is less arbitrary in this aspect, because one can start 
from earlier stage of the evolution, in which the choice of uniform 
distribution in P 2 /P 1 and log 10 ai is better. 


8 CONCLUSIONS 


We have studied the migration of two super-Earths embedded in a 
protoplanetary disc. The forces acting on each planet are given by 


analytical prescriptions (Paardekooper et al. 2011 

Tanaka & Ward 

2004,|Fendyke & Nelson 2014, Dittkrist et al. 201 

4). It enabled us 


to study the orbital evolution of many initial configurations over the 
whole disc life-time. We showed that there are regions in the disc 
in which the migration is divergent. One of the regions corresponds 
to the outward migration zone, which is placed in the shadowed 
part of the disc (radii between ~ 1 and a few astronomical units) 
as well as in a transition zone between the irradiated and the shad¬ 
owed parts of the disc. The second region of the divergent migra¬ 
tion (at around of 0.5 AU) is a result of using more realistic opacity 
law (S emenov et al. 2003) instead of a simpler law (Ruden & Pol-| 
|lack|1991| >. The third zone of the divergent migration (r < 0.1 AU) 
corresponds to a region in which the opacity is dominated by gas. 
Other regions in the disc corresponds to the convergent migration. 
In general, during the evolution, a given system spends some part 
of the time in regions of the convergent migration (the period ratio 
decreases and can be halted at one of the resonant values, i.e., 2 / 1 , 
3/2, 4/3, etc.) and some other part of the time the system spends in 
regions of the divergent migration (the period ratio increases). The 
final configurations depends on the initial orbits as well as on the 
masses of the planets. 

We performed 3500 simulations for planets’ masses and ini¬ 
tial orbits chosen randomly. Almost all the systems ended up as 
periodic configurations. We found that the period ratio axis can 
be divided into several parts. The first part (P 2 /P 1 > 2 . 12 ) corre¬ 
sponds to non-resonant configurations, although one of two reso¬ 
nant angles of 2:1 MMR oscillates. Systems located in the second 
part of the period ratio axis (P 2 /P 1 E [~ 1.72, ~ 2 . 12 ]) are involved 
in 2:1 MMR. Both resonant angles librate and the amplitudes of 
the librations depend on P 2 /P 1 . For P 2 /P 1 ~ 2 the amplitudes are 
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« 0. The amplitudes of both angles are monotonically increasing 
functions of P2/P1 when P2/P1 >2 and monotonically decreasing 
functions of the period ratio when P2/P1 < 2. For P2/P1 ~ 1.72 
the amplitudes of both angles reach 360 degrees, which means that 
the system is located at the border of 2:1 MMR. For P2/P1 ~ 2.12, 
the amplitude of the second angle reaches 360 degrees, while the 
amplitude of the first angle is still lower that 10 degrees. The third 
part of the period ratios axis corresponds to 3:2 MMR. The range 
of period ratios of this resonance is [~ 1.4, ~ 1.7]. Similarly to 
2:1 MMR the amplitudes of the resonant angles librations are « 0 
for P2/P1 ~ 1.5 and they increase when P2/P1 deviates from 1.5 
in either of the sides. At the borders of this resonance, the am¬ 
plitudes reach 360 degrees. Further parts of the period ratio axis 
corresponds to 4:3 and 5:4 MMRs. The amplitudes of the resonant 
angles librations have similar functional form as for the two previ¬ 
ously discussed resonances. 

The general conclusion is that almost all the systems with 
P 2 /P 1 ^ 2.12 are involved in one of the first order MMRs. Never¬ 
theless, for systems whose period ratios are far from (/?+!)//?, the 
ranges of e cos 05 and e sin05, in which the resonant angles librate, 
are very small (~ 10 -4 -i- ~ 10 -3 ). It means that even small pertur¬ 
bations of the orbits (e.g., due to turbulences in the disc or inter¬ 
actions with planetesimals) would probably change the behaviour 
of the resonant angles from librations to rotations and also move 
the system out of the periodic configuration. The fact that most of 
the systems which underwent the smooth migration (convergent or 
divergent) evolve periodically might be tested with observations. If 
a given observed system was found to be close to a periodic orbit, 
it would mean that the system unlikely experienced strong pertur¬ 
bations at early stages of the evolution. This kind of studies could 
give the insight into the magnitude of turbulences in the disc as well 
as number and sizes of planetesimals. 
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